samplesplit=function(obj,col,num) {
  library(matlab)
  result=data.frame(matrix(,nrow=num,ncol=3))
  colnames(result)=c("bin","coef","cse")
  
  for (i in 1:num) {
    result[i,1]=LETTERS[i]
    
    if (i==1) { result[i,2]=obj$coefficients[which(rownames(obj$coefficients)==paste(col,LETTERS[i],":","B",sep=""))]
    } else {
      result[i,2]=obj$coefficients[which(rownames(obj$coefficients)==paste(col,LETTERS[i],":","B",sep=""))]-
        obj$coefficients[which(rownames(obj$coefficients)==paste(col,LETTERS[i],":","W",sep=""))]
    } 
  } 
  
  for (i in 1:num) {
    if (i==1) { result[i,3]=obj$cse[which(names(obj$cse)==paste(col,"A:B",sep=""))] 
    } else {
      pos1=which(rownames(obj$coefficients)==paste(col,LETTERS[i],":","B",sep=""))
      neg1=which(rownames(obj$coefficients)==paste(col,LETTERS[i],":","W",sep=""))
      test=rep(0,nrow(obj$coefficients))
      test[pos1]=1
      test[neg1]=-1
      result[i,3]=sqrt(test%*%obj$clustervcv%*%test)
    }
  } 
  result$lb=result$coef-1.96*result$cse
  result$ub=result$coef+1.96*result$cse
  result
} 


samplesplit.bh=function(obj,col,num) {
  library(matlab)
  # implicit: reference level is A:W
  result=data.frame(matrix(,nrow=num,ncol=3))
  colnames(result)=c("bin","coef","cse")
  for (i in 1:num) {
    result[i,1]=LETTERS[i]
    
    if (i==1) { result[i,2]=obj$coefficients[which(rownames(obj$coefficients)==paste(col,LETTERS[i],":","BH",sep=""))]
    } else {
      result[i,2]=obj$coefficients[which(rownames(obj$coefficients)==paste(col,LETTERS[i],":","BH",sep=""))]-
        obj$coefficients[which(rownames(obj$coefficients)==paste(col,LETTERS[i],":","W",sep=""))]
    } 
  } 
  
  for (i in 1:num) { 
    if (i==1) { result[i,3]=obj$cse[which(names(obj$cse)==paste(col,"A:BH",sep=""))] 
    } else {
      pos1=which(rownames(obj$coefficients)==paste(col,LETTERS[i],":","BH",sep=""))
      neg1=which(rownames(obj$coefficients)==paste(col,LETTERS[i],":","W",sep=""))
      test=rep(0,nrow(obj$coefficients))
      test[pos1]=1
      test[neg1]=-1
      result[i,3]=sqrt(test%*%obj$clustervcv%*%test)
    }
    
    
  } 
  result$lb=result$coef-1.96*result$cse
  result$ub=result$coef+1.96*result$cse
  result
}


barplot.arrows=function(df,num,ylims,names,main,ylab,clr) {
  barplot(df$coef,width=1,space=.25,ylim=ylims,col=clr,
          names.arg=names,
          main=main,cex.main=1.5,cex.axis = 1.2,cex.names = 1.1,
          ylab=ylab)
  xpos=1:num*1.25-1/2
  par(new=T)
  arrows(xpos, df$lb, xpos, df$ub,length=.05,angle=90,
         lwd=2,col="black",code=3)
}